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We solve the Bethe-Salpeter equation (BSE) for a system of a heavy quark-antiquark pair inter- 
acting with a Poincare invariant generalization of screened linear confining potential. In order to 
get reliable description the Lorentz scalar confining interaction is complemented by the effective one 
gluon exchange. Within presented model we reasonably reproduce all known radial excitations of 
the vector charmonia. We have found that J/SP is the only charmonium left bellow naive quark- 
antiquark threshold 2m c , while the all excited states are situated above this threshold. We develop 
a method which is enable to provide solution of full four dimensional BSE for the all excited states. 
We discuss the consequences of the use of the free propagators for calculation of excited states above 
the threshold. The Bethe-Salpeter string breaking scale /i ~ 350MeV appears to be relatively larger 
04 ' then the one defined in various potential models fi ~ !50MeV. 
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I. INTRODUCTION 



Excited meson spectroscopy is a keystone experimental output, which is essential for understanding of quark- 
antiquark interaction. The dynamics of meson constituents is driven by solely known strongly interacting quantum 
field theory- quantum chromodynamics. The explanation of quarks and gluons confinement is one of the great 
i-C ' challenge of theory of strong interaction. The same confining forces are responsible for a large degeneracy which 
emerges from the spectra of the angularly and radially excited resonances. As reported in [H-Q such degeneracy is 
pi | observed in pp annihilation by the Crystal Ball collaboration at LEAR in BERN [5] . Similar can be deduced from 
heavy quarkonia production in e + e~ angulations. Recently, BABAR, Belle, BESS and LHCb experiments continue 

• collection of various meson experimental data. 

Following ideas of Ref. @ , confinement in heavy flavor sector has typically been associated with a linearly rising 
£SJ ■ potential between constituents Q • The spin degeneracy observed in heavy quarkonia spectrum tell us that the main 
^ | part of confining interaction should be largely spin independent. A various lattice fits leaded to various predictions 
for the potential between quark-antiquark states. A well known is Cornell parameterization of the static Wilson loop 

• derived potential @, Q, 



V(r) = -a/r + ar (1) 



which appeared to be suited for description of the first few excited mesons. 
Recently in Refs. [IM1 

it has been found that the meson spectroscopy is better described by " confining" potential 
which is bounded from above. While in the absence of dynamical quarks the nonrelativistic static potential can 
eventually show up a linear asymptotic, it should be flattened in large distance due to the string breaking associated 
with light hadron productions. In this respect, the exponential potential can be regarded as a screened version of the 
linear potential. The screening effect should be universal, e.g. scheme and gauge independent property of low energy 
QCD as the creation of the light quark-antiquark pairs is energetically favorable and pions and other light mesons 
is observable fact. The screening could be included in order to explain observed hierarchy of radially excited heavy 
mesons. Obviously, the spectrum deviates from the linear Regge trajectory, for a proposal of 5s and 6s charmonium 
candidate see [l3j], however recall, the deviation from the linear Regge trajectory is expected in the light meson sector 
as well Q. 

The Coulombic part of the potential is naturally expected in QCD, however as shown recently, it can be suppressed 
in charmed meson sector by not yet understood mechanism [l5j . Furthermore, it is related with gluon exchange then 
screening of Coulomb potential is expected as well. This a bit more subtle matter could be related with soft gluon 
mass generation [16l |17| through the Yang-Mills Schwinger mechanism. Actually the dynamical gluon mass generation 
is suggested by finiteness of the lattice gluon propagator in the deep infrared, for the recent lattice data on gluon 
propagator in Landau gauge see [3]. We expect, the screening mass characterizing string breaking and the soft gluon 
mass have similar size ~ A-qcd- 
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For bound states which lie above the naive quark-antiquark threshold, the BSE kernel becomes singular, which 
causes that many usual numerical treatments fails and/or become impossible in practice. To authors knowledge, 
plethora of solutions have been obtained in relativistic quantum mechanic [l^] or by solving various phenomenological 
3D reduction of BSE [2(| [Hj], however comparison with the original BSE is completely missing. The way of three- 
dimensional reduction of the BSE is rather arbitrary and difficult to control without keeping the solution of the original 
BSE. Actually, choice of 3d equation has a large effect already for the ground states. In the study 20] 200MeV 
difference between Equal Time (ET) and certain Quasipotential Equation has been found for the charmonium ground 
state I = 0. This increases slightly for higher I as well. The differences between Schroedinger and 3d approximated 
BSE spectra appears to be more significant for higher excited states, where following ref. |2(| the difference has 
been estimated Mqm — Met — 300 MeV for I = 4, 5 charmonia. The comparison between various relativistic 
approximations has not been studied for radially excited states, however we expect very similar effect there. From 
all this is apparent that relativistic covariance is important for charmonium and it will be useful to have the solution 
without any 3d approximation. Due to this reasoning, using single component approximation, but keeping the full 
four dimensionality of BSE we study the effect of retardation in the slopes of radially excited vector mesons. The 
experimental knowledge of heavy vector quarkonia is the main reason why we concern spin 1 mesons as a first (as 
time is being the pseudoscalars have been computed already in [2l[ ) . 



II. KERNEL FOR BETHE-SALPETER EQUATION 



In the Quantum Field Theory the two body bound state is described by the three-point bound state vertex function, 
or equivalently by the BS amplitude. Both of them are solutions of the corresponding covariant four-dimensional BSE 
[23j . In principle, common framework of Schwinger- Dyson and BSEs offers unique Poincare invariant generalization of 
quantum mechanical picture sketched above, however a practical solutions are always incomplete due to the truncation 
of the equations system. Due to this fact we rather phenomenologically estimate what should be the form of the BSE 
kernels here. In this paper we use " confining" interaction kernel of the form 

V a (q) = f 2 ° 2 , 2 ; (2) 
{q 2 - M ) 

which is certain screened form of ~ I/9 4 scalar interaction. The remaining considered part of the interaction kernel 
has the Dirac decomposition identical with the one gluon exchange ~ 7q^K7„.^,- 
For clarity we write down the BSE completely here 



S- 1 ( q + P/2)xi P ,q)S-\q + P/2) = -i J -^^ x (k,P) lv G^{k-q) 



(2i 



X (k,P)V s (k-q): 



G^[k - q) = g»"V v = - g2 ^T 2 5 (3) 
(k - q) 2 - fi 2 g 

Two scalar functions V s and V v in Q completes our simple Poincare invariant generalization of quantum mechanical 
potentials. Here, clearly G^ u represents effective gluon propagator in Feynman like gauge, where the effective soft 
gluon mass \x g has been introduced. The double pole scalar interaction V s leads to regular exponential potential in 
the position space. Actually, in heavy quark limit one can consider three dimensional potential 

V^(k) = [°° dr ^Mkr) v{r) 







(4) 
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where the potential in position space reads 



(k 2 +^ 



V(r) = -a . (5) 
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Thus in certain sense the BSE model considered in this paper represents relativistic generalization of the models 
considered in (ToL [Tl|. S in Eq. stands for charm quark propagator, which in our simplest approximation is taken 



as 



S =/ — m c , 

and xv represents the Bethe-Salpeter wave function which has the general form: 

Xv(q,P) =4xvo+ /Pe.qxvi + fc-qXV2 + e.qxv3 

+ [4, /P]XV4 + [4, ff\XVh + [A, f]XV6 + H5 fXV7 , 



(6) 



(7) 



with t„ = e^q u P a e^ e 2 = -1, e.P = 0. 

Munczek and Jain _24j have shown that V0 component is dominant for the all ground state mesons, which dominance 
is particular for mesons made form heavy flavor (anti)quarks. The same is valid for the case of all pseudoscalar radial 
excitations [29j , which strongly argue for that if there is a dominant component it stay to be dominant one for higher 
excited states independently on the meson spin. Therefore we assume the same apply for excited vectors here, and 
we neglect all other components in presented study. 

Within the approximation the Bethe-Salpeter equation in the rest frame reads 



Xvo(pe,P) 



P 2 /4 



where 



-Pi 



d A k E 
(2^ 



Tfl 



M 2 /4) 2 + ql 



(8) 



-2V v + V s } X vo(k E ,P) ; 



V v 



v s = 



q% + [i 2 



C 



(9) 



where we have performed Wick rotation into the Euclidean space. For Euclidean momenta we use the convention 
k E = {ki, k), k E = k\ + k 2 , while the total momentum is kept timelike P 2 — —P E — M 2 as required for the bound 
states. An extra problem arises for bound states which are heavier then sum of constituents quark masses. As we are 
using propagators with single real poles, the BS wave function becomes singular since it is proportional to the product 
of the quark propagators. For the vertex function, the threshold like singularity should appear for the solution with 
P 2 > 4m 2 . Due to this fact the numerical integration requires special numerical care and we found very advantageous 
to define the following auxiliary function 



A{ PEl P) = X vo(PE,P)G A (k,P) 



G A (k,P) = 



Pi 



P 2 /4 



-Pi 



m 2 + M 2 /4) 2 + p 2 4 M 2 + ■ 



(10) 



where e is a small regulator mass satisfying e 4 << p 2 M 2 . Within the numerics the regulator can be limited to zero. 
In limiting case we should add the residuum contribution stemming from the propagator pole. In the presented study 
we keep the regulator always finite and we neglect the residuum as well. It avoids complicated analytical continuation 
to the Minkowski space in this case. 

Finally we integrate over the angles of 3d momentum subspace. Explicitly written the BSE (|8|) reads 



■^r ^ dkk 2 A(k E , P)G A (k,P) J j^t [-2V V + V s ] 

which can be finally rewritten as 

dk A / dkA(k,P)G A (k,P) [K s + K v ] ; 

-oo JO 



A(p,P) 



(11) 



(12) 



C 



(2tt) 3 [ k % + p %- 2k 4P4 + /i 2 ] - 4k 2 p 2 



K v = - 



g 2 k 



(2tt) 3 p 



-In 



k 2 + p 2 - 2k 4 p 4 - 2kp + \i 2 



k 2 +p 2 - 2fc 4 P4 + 2kp + p 2 
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FIG. 1: Comparison of BSE solution with PDG data for vector charmonium. 



where the functions K,, 



stem from Q integration: 
(2tt)4 



Vv = K v 



■.iii 



(2tt)4 



V s = K s 



(13) 



III. NUMERICAL SOLUTION OF THE CHARMONIUM BSE 



It is well established that in unconfining theory a bound states spectra obtained through the BSE and through the 
corresponding Schroedinger equation mutually agree up to the small relativistic correction. In opposite, the solution 
of BSE for excited states, which lie above constituents particle threshold represents rather difficult numerical problem 
and no comparison is known in the literature, at least to the author. 

Due to the presence of the kernel singularity more or less standard matrix methods [25| fail since the inversion of 
numerical matrices is not possible. Also we do not explore more or less conventional expansion into the orthogonal 
polynomials which loses its efficiency when, as one expects, relatively large number of polynoms is necessary. Instead 
of, we rather solve the full two dimensional integral equation by the method of simple iterations. For this purpose we 
discretize P 2 and step by step we are looking for the solution of the BSE with given Pf . 

The BSE for bound states is a homogeneous integral equation and it satisfies usual normalization condition. Instead 
of using this, to achieve a good numerical stability of the iteration process we implement normalization condition 
through an auxiliary function X(P) and solve the following equation: 



A(p,P) = X(P) J dk A J dkA(k,P)G A (k,P)K E . 
The following has been found as a particular useful choice for the function X(P) 

A" 1 (F)= J dki J dkA(k,P) 2 f(k,P), 



(14) 



(15) 



where arbitrary positive weight function / was chosen to be Gaussian in k and k^. Implementation of such X(P) makes 
BSE nonlinear but mainly numerically stable. Clearly the BSE solution has been identified when X(P) = 1 and when 
the difference between consecutive iterations vanishes at the same time. We found that these two conditions happen 
simultaneously, while for other values of parameters P, X(P) =/= 1 the numerics do not provide vanishing difference 
between iterations. 

Numerical convergence of A — > 1 depends on the density of the integration points in the important domain of 
momenta. Recall here that even for 1-dim reduced BSE see 27| relatively large number of integration point was 
required in order to get precise solution. In presented study we are dealing with 2-dim integral equation with 
principal value integration, thus to achieve the same acuracy it necessarily enlarge number of integration points. 
Compromising between computational time and an estimated numerical error we are satisfied with the use of maximum 
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M(32) 


A(32) 


cr(32) 


M(58) 


A(58) 


cr(58) 


M(72) 


A(72) 


«r(72) 


M(96) 


A(96) 


o-(96) 


3.717 


0.995 


1.2E-06 


3.642 


1.398 


0.006 


3.639 


1.444 


0.008 


3.636 


1.458 


0.008 


3.832 


1.015 


1.5E-05 


3.746 


1.174 


0.001 


3.745 


1.200 


0.002 


3.733 


1.223 


0.002 


3.914 


0.994 


2.0E-06 


3.865 


0.998 


2.5E-07 


3.867 


1.000 


1.7E-09 


3.855 


1.017 


1.8E-05 


4.040 


0.993 


2.5E-06 


3.999 


1.027 


4.5E-05 


3.996 


0.979 


2.6E-05 


3.980 


0.985 


1.23E-05 


4.151 


0.983 


1.7E-05 


4.061 


1.020 


2.5E-05 


4.049 


0.983 


1.7E-05 


4.033 


1.01 


2.24E-05 


4.289 


1.012 


1.0E-05 


4.177 


1.024 


3.7E-05 


4.168 


1.002 


3.5E-07 


4.149 


1.00 


5.24E-06 


4.458 


1.024 


3.7E-05 


4.277 


1.030 


5.7E-05 


4.262 


1.002 


3.5E-07 


4.243 


1.01 


2.03E-05 


4.604 


1.025 


4.0E-05 


4.402 


1.013 


1.1E-05 


4.340 


0.977 


3.1E-05 


4.365 


1.02 


1.3E-05 


4.840 


1.017 


1.8E-05 


4.549 


0.992 


3.5E-06 


4.530 


0.980 


2.4E-05 


4.503 


1.02 


1.4E-05 


5.008 


1.021 


1.8E-05 


4.690 


0.998 


1.4E-07 


4.671 


1.007 


3.98E-06 


4650 


1.005 


5.1E-06 


5.318 


0.989 


7.00E-06 


4.890 


0.99 


1.4e-07 


4871 


1.01 


2.07E-05 








5.518 


1.009 


5.5E-06 


5.055 


0.998 


1.4E-07 


5.030 


1.009 


5.09E-06 








5.934 


1.018 


2.1E-05 


5.337 


1.013 


1.0E-05 


5.308 


1.01 


1.20E-05 








6.171 


1.015 


1.5E-05 


5.527 


0.987 


9.8E-06 


5.496 


1.02 


2.60E-05 








6.728 


0.998 


1.6E-07 


5.909 


0.97 


3.5E-05 




















6.140 


1.00 


2.5E-06 















TABLE I: Spectrum obtained for various number of mesh points N = 32, 58, 72, 96 used at each single integral. The density of 
mesh point is regulated by A = lOOOGeV and infrared regulatore = 0.03GeV^ is used, A and the iteration error a are displayed 
for completeness. There is one more state observed at 3.880 (0.998, 2.5i? — 07) for N=58, which is skipped in the table for 
better comparison. Mj<s, = 3098 is used as a fit. 



Nki * -/Vk = 184 * 96 Gaussian points. In order to ensure the numerical stability we varied integration volume (in 
p-space) and the number of points as well. In order to get optimal density of integration points the upper boundary 
A is introduced in momentum integrations. More precisely, for the variable k we map single integration interval (0, 1) 
into (0,A) by simple rescaling (and (—1,1) into (—A, A) for variable k^). To assure convergence A —¥ 1 the density 
is varied such that for lower states smaller A is used, while highly excited states are better obtained within higher 
A, with some cost of numerical precision in the later case. The numerical dependence on N is shown in the Tab. 1. 
We also compare with the experimental data in the Tab. 2 and Fig. 1. For purpose of brevity we do not discussed 
all numerics here leaving more detailed discussion elsewhere. However, we can mention that as a numerical test we 
have checked the numerical method against the scalar models [2q . where the resulting BSE spectra have been already 
obtained by different (and in fact more accurate) method. Thus we solve 2dim BSE on a grid given by the component 
of relative momentum k^ and k within the procedure described above, which provides the spectra within ~ 1/N% 
accuracy (N is number of points in the single, say k- integration). As the regulator e is implemented here, we expect 
the accuracy of excited vectors spectrum studied here could be comparable with the one of scalar models solutions. 
In our case we estimate the 1 % numerical error in determination of bound state masses. Such numerical error is 
considerably smaller then the expected amount of energy shift due to the open threshold effects ( effect associated 
with D meson production, not be confused with unphysical quark-antiquark threshold. An open charm is impossible 
to incorporate in the formalism of homogeneous BSE and it is completely ignored here). 

For clarity we restrict ourself to the dynamics of charmonia in this paper. The numerical values of the models are 
the following 

C* = 5A18GeV 2 ;a s = g 2 /{An) = 0.2 (16) 
m c = 1.5615GeV; fi = \i g = 364.35MeV 

and we use the experimental J/\f r mass to fine tune the correct scale at the end (in units where m c = 1.5 ; C — 5 
we have M(J/\&) = 2975). In order to improve the numerical stability the infrared regulator in (|10p was adjusted as 
e = O.QSGeV (for A = lOOOGeV^, e value is taken slightly larger then the smallest interval between integration points 
of the variable kn). 

As one see from the Table I, numerically A ^ 1 for the first two states. In fact A do not cross unit value, instead 
it posses positive local minimum which is stated in the Table. Further, when comparing with the experimental data, 
there are two more BSE solutions with masses in between ip and ip , while the rest of calculated excited states quite 
nicely agree with the data. We have not found simple way (by varying parameters C, a, mu, m c ) to exclude these two 
additional states from the solutions and we suggest the reason why they are here in the Section bellow. 



6 



M m 


PDG 


n, I 


3097 


3097 


Is 


3640 


3686 


2s 


3730 


3772 


Id 


3850 






3980 






4030 


4039 


3s 


4150 


4153 


2d 


4240 


4263 


4s 


4360 


4361 


3d 


4450 


4421 


5s 



TABLE II: Comparison with PDG data (second column) and calculated spectrum. Quantum numbers correspond with assumed 
quantum mechanical assignment [12]. 



IV. CONCLUSION AND DISCUSSION OF THE RESULTS 



We have formulated Lorentz covariant model for vector quarkonia, which is based on the BSE with phenomenological 
kernel. With aforementioned exception of two additional states, the resulting spectrum is comparable with the 
experiments whenever the experimental data are available. The agreement between our results for higher states and 
the one measured in the experiments is impressive, remaining difference between theory and experiments is due to 
the approximations, e.g. due to the interplay of quark and D, D* mesonic degrees of freedom, such couple channel 
effects are difficult to incorporate into the Bethc-Salpeter analysis presented here. 

The ground state is situated near bellow the naive quark-antiquark threshold, while the all excited states lie 
above this threshold. We argue two more states which appear are the artefact of inappropriate usage of free quark 
propagators. The question of confinement is beyond scope of the presented work, however we expect some changes 
when confinement is correctly incorporated. First of all, we expect the quark propagators should not have a free 
particle pole and therefore the BSE could not posses ordinary threshold singularity in this case. In the paper we are 
dealing with BSE where the propagators describe free -instead of confined- quarks. Therefore the threshold singularity 
unavoidably appears as an artefact here. The lowest lying excited charmonia are the one closest to the naive quark- 
antiquark threshold and we naturally must expect some defect in the calculated spectrum. Here, very likely it leads 
to mentioned appearance of two more excited states that we cannot find in the PDG, (for a recent attempt to find 
more definite answer see [29|). 

The model presented is very simple: The BSE kernel consists from vectorial effective one gluon exchange and from 
the scalar infrared enhanced -double pole- interaction. A possible vector-scalar admixture of "confining" interaction 
has not been considered and we expect it must be small in order to suppress large hyperfine splitting. The string 
breaking mechanism is incorporated through the screening mass fi which is found to be comparable to Aqcd in our 
BSE study. Presence of Lorentz scalar interaction complemented by vector-vector interaction kernel appears to be 
very important part of the model. No individual -scalar nor vector -interactions provide quarkonium spectrum [28j . 

Nonrelativistic quantum mechanical limit of the scalar interaction is given by the exponential potential. Conse- 
quently the spectrum of radially excited states do not correspond with linear Regge trajectory but the gap between 
the states increase with the mass of the bound states. This is in very good agreement with recent experiments. The 
intercept -the gap between J /iff and the first excited state *ff(2S) is driven by the interplay of the strength of vectorial 
and scalar interaction. As we have checked numerically, coulombic term can be regarded as a perturbation only for 
existing energy levels for which it slightly shifts existing energy levels. On the other side, its omission would lead to 
the appearance of new states bellow existing 2S mass. We found that the correct adjustment of the energy intercept 
is major effect of the effective one gluon exchange interaction. The running coupling is fixed here and its numerical 
value differs significantly from the one known from quantum mechanical phenomenology (it is more then two times 
smaller). At present stage we have no simple explanation of this fact. If our results are considered seriously then 
it can suggest that average of typical square of gluon four-momenta inside quarkonium can be considerably larger 
then naively expected from Schroedinger equation. A possible explanation is that the soft gluon mass is larger then 
the value we consider here , fi g = 700 — 800MeV, however QM mechanical limit has not yet been studied in such 
case. Further, comparing with the expected quantum mechanical limit then Bethe-Salpeter strin g b reaking scale 
fi ~ ShOMeV appears to be relatively larger then the one used in potential model /i ~ 150MeV [T51 1 111 ], 
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There is another challenging aspect of this problem, as most of the SDE-BSE studies rely on the ladder truncation 
of the equations system. After the inclusion of running quark masses, renormalization wave functions, the techniques 
can be very useful in calculation of heavy-light and light flavored meson as well. We also expect that the knowledge 
of off-shell behavior (e.g. q.P dependence) of the amplitudes is important in various hadronic processes. Due to this 
it is worthwhile to extend our study to the more complete calculations, which open possible first principle calculation 
of production and decay mechanisms of heavy flawour hadrons. 
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